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PREFACE 


This report was prepared by A.R.T. Research 
Corporation, Los Angeles, California, under 
Contract NAS3-14400 and was funded by the 
National Aeronautics and Space Administration- 
Lewis Research Center, Cleveland, Ohio. Inclusive 
dates of research were 24 June 1970 through 
4 December 1970. The NASA Project Manager for 
this work was Mr. Millard L. Wohl. 

This report comprises two (2) volumes; Volume I- 
Summary Report covers the theoretical basis for 
the FASTER-III computer program and results for 
sample problems; Volume II - Users Manual gives 
detailed operational instructions for the computer 
program . 
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ABSTRACT 


This volume outlines the theory used in FASTER-III, 
a Monte Carlo computer program for the transport 
of neutrons and gamma rays in complex geometries. 
The program includes the treatment of geometric 
regions bounded by quadratic and quadric surfaces 
with multiple radiation sources which have a speci- 
fied space, angle, and energy dependence. The 
program calculates, using importance sampling, 
the resulting number and energy fluxes at specified 
point, surface, and volume detectors. It has 
the additional capability of calculating the 
minimum weight shield configuration which will 
meet a specified dose rate constraint. 

Results are presented for sample problems involving 
primary neutron and both primary and secondary 
photon transport in a spherical reactor-shield 
configuration. These results include the optimiza- 
tion of the shield configuration. 

The users manual for the FASTER-III program is 
contained in a companion volume. 
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Section 1 


INTRODUCTION AND SUMMARY 

The origins I PASTER program. Reference 1, contained a number 
of new techniques which provided the capability of obtaining 
accurate radiation levels at specified points in complex 
geometries. The use of this program by NASA-Lewis Research 
Center and other Government facilities and contractors in- 
dicated the need to broaden the overall program capabilities, 
automate the importance sampling, increase the computational 
efficiency, and revise the users manual. This revised program 
has been designated FASTER-III to distinguish it from earlier 
versions . 


A specific program capability developed for NASA-LeRC permits 
a calculation of minimum weight shield configurations for 
mobile nuclear reactor applications, e.g., nuclear propulsion 
for aircraft, surface effect vehicles, and space craft. 

The basic Monte Carlo transport method was extended to include 
a calculation of partial derivatives of the radiation fluxes 
with respect to specified shield dimensions. These derivatives 
are then used to define exponential relationships used in 
the shield optimization procedure. This optional program 
feature is described more completely in Section 2. 

A number of program revisions had also been made by A.R.T. 
Research Corporation for various customers and to provide 
an internal capability for solving a variety of radiation 
transport problems. These revisions are included in the 
FASTER-III program. Particularly noteworthy are the following: 

(l) A calculation of optimal importance sampling 

parameters based on partial derivatives of the 
variance (Section 2.3). 
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(2) The acceptance of data in either fixed or 
variable field formats including the ANISN-DTF 
format for neutron cross sections. 

(3) The calculation of time-dependent neutron 
and photon transport (using time . moments and/ 
or time intervals) including an optional 
exponential atmosphere. 

(4) The improvement and addition of importance 
sampling models with the various importance 
sampling parameters built into the program. 

Various program features are described in References 2-6. 

The application of the FASTER-III program to a shield optimi- 
zation problem is discussed in Section 3. The problem in- 
volved a spherical reactor-shield configuration and included 
primary neutrons and both primary and secondary photons. 
Conclusions and recommendations are presented in Section 4. 

Volume II (Users Manual) presents the detailed description 
of the FASTER-III program along with all the instructions , 
for operation on the IBM 7094, UNIVAC 1108, CDC 6600, and 
IBM 36O-0S (single or double precision) computers. 


2 



Section 2 


ANALYSIS 

The techniques used in calculating optimum shield configurations 
and optimum importance sampling parameters are summarized below. 
The discussion is given in three parts: derivatives of fluxes 

with respect to shield dimensions, optimization techniques, 
and derivatives of variance with respect to importance sampling 
parameters. The basic Monte Carlo techniques assumed in this 
discussion are summarized in Appendix A. 

2 . 1 Shield Dimension Derivatives 

The dose rate at a point 'detector jr for a specified reactor 
shield configuration is written as: 

J 

D (z) = (U 

j = L 

where J is the total number of energy groups for both neutrons 

and photons (including secondaries), is the particle flux 

in the jth energy group, and R. is the conversion factor from 

J 

flux to dose rate. The rate of change of the dose rate with 
respect to a shield dimension is simply 


dD (y) _ 
St i 



£ = 1, 2, ..., L (2) 


where L is the total number of shield dimensions and t^ is the 
value of the ^th shield dimension. 
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The equation used by the program for determining the flux is 
written as : 


0j(z)- n 


N 


EE 

n= l k 


S jkn^-kn.) V^kn- 




l^tal 


( 3 ) 


where N is the total number of histories tracked via the Monte 
Carlo method, k is the number of particle collisions, 

-kn is the P° sition of the k th collision of the n th history, 
S jkn^-kn ) the number of particles in the jth energy group 
emerging from in the direction u^ of the detector per 
unit solid angle, and ^(z^,^;) represents the material and 
geometric attenuation kernel for particles in the jth energy 
group going from z^ to the detector. 


The partial derivative of the flux with respect to a shield 
dimension is simply: 


(z) 


1 

N 






The summations are a minor part of the calculation. Therefore, 
the notation is simplified by concentrating on the elements 
in the summation 

“Tt J~ = ~btj [ S jkn^-kn ^ K j^-kn^^] ^ 

where 9^^ represents the contribution to the flux in the jth 
energy group from the k th collision of the nth history. 
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This equation is rewritten as 


d 9 ,ikn _ Q 


atg jkn dtg 


In 


^jkn^— kn^ kn'^ 


( 6 ) 


0 jkn 


% ln S lkn^ to ) + -fy ln K j ( ^kn^ ) 


The second term in brackets involves the attenuation kernel 


exp 


V 2 W£> = 


M 


m=l 


S a . 
m jm 


( 7 ) 


where M is the total number of regions traversed from z to 

~k n 

the detector, s m is the path length for the m th region traversed, 
a . is the total cross section of this region for particles 
in the jth energy group, and s is the total distance from 
to the detector, i.e.. 


s = 


M 


m=l 


(8) 


m 


A substitution of this kernel gives: 

M 


sir ln s 


3t i 


M 


^Lrf S m °jm ~ 2 In ^^ s m 


m=l 


m=l 


M 


E 

m=l 


ds 


jm 


m 2 ^ w ~m 


ds r 

rn= L bt£ 


l 


s 


m=L 


m 


M 

= " £ + 


2 N as m 


at 


l 


(9) 
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The partial derivative of the partial path length s m with 
respect to the shield dimension t^ is zero unless the mth 
region traversed is affected by a change in t ^ • In particular, 
if t| is a characteristic dimension of the region, i.e., 
its thickness, then 


as 


m 




knm 


knm 


— kn*-knm 


( 10 ) 


where is the cosine of the angle measured from the surface 

normal > with which the particle crosses the boundary of 

the region. 


In the strict sense, the change of the dimension of one shield 
region can affect other shield regions. In particular, for 
a spherically symmetric reactor-shield configuration, an increase 
in the thickness of a shield region forces a movement of all 
shield regions having a larger radius. The inclusion of these 
effects in the above equation unnecessarily complicates the 
analysis and the calculations. The primary effect of changing 
a shield region dimension is to change the number of mean free 
paths which particles have to traverse in reaching the detector. 
Therefore, in calculating the derivatives, only the effect of 
the material attenuation is treated. 


The derivatives at a specific boundary crossing m' then simplify 
to: 

■ M 


r- l” VSto-z) = + §> 


as 


m 


7 


at, 


m-l 


1 


/ 2 \ 1 / r\ , 2 \ 1 

= - ("jm- + I>K - <° + 


knm' 


s 7 -M-. 


knm 1 


= " ‘jn/^knm' < U > 

where m' is the index of a region having t* as a dimension. 
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The partial derivatives of the particle weight with respect 
to the shield dimensions -- the first term in brackets in 
equation 6 -- are zero at the point of origin of all primary 
particles. For subsequent particle collisions, the derivatives 
are calculated using the relationship between particle weights 
on subsequent collisions: 


Sjkn^— kn ^ 


^ °i,k-l,n^kn ^ K i^- 2 k-l,n , ~kn ^ T i j kn*— kn*-kn ^ 


p kn ^ -kn ^ 


-kn 


~kn ~k-l,n 
2 kn -k-l,n 


( 12 ) 


where S* n (x^ ) is the number of particles coming out of 
the previous collision point in the direction v^ and in the 


i th energy group, K i (z^_ 1 n ,z kn ) is the attenuation kernel be- 
tween particle collision points, T^j (z^ ,v^ n * u^ ) 


is the scatter- 
ing kernel for transfer of particles from group i to group j, 

* / \ 

and is th e probability density function used in select- 

ing the collision point. 


A straightforward substitution gives 




In 


S S i , k - 1 , n ^ -kn ^ K i ( -k - 1 , n ' -kxK) 1 ^kn "4mHhn ) 


p kn ^ -kn ^ 


( 13 ) 
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After some manipulation, this reduces to 


% ln s jkn ( H to ) = TS~ " E V ijkn 

* jkn l -kn' i 


SU ln S l,k-l,n<2ta> 


+. In K. (z. . , z, ) - — ^ — In r*(z, ) 

i — k-l^n^-kn ' $ta Tarf-kn' 


at 


t 


W 


where 


V 


= S i>k-l,n^-kn^ ^^-k-l.n^-kn ^ T i.j ^kn'-kn* -kn l 


ijkn 


Pkn^ta,) 


(15) 


The first term in brackets in equation 14 is the same partial 
derivative for collision k-1 as the partial derivative now 
being calculated for collision k. Therefore, it is known, 
either identically zero for k=0, or as determined from equation 14 
for k>0. The second term in brackets in equation 14 is similar 
to the second term in brackets in equation 6 and is therefore 
determined by equation 11. The last term in brackets involves 
the definition of the probability density function used to 
select the collision point z^ . 

The probability density function for a collision point has the 
form 


vL(^) - OO A(s)a(s) exp 

. -J 

r 

o 

a(s ' )ds * 


A(s ' )a(s ',) exp 

[i 

^* s a(s")ds 

3 ' 

1! 

ds * 
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’K’ 

where i s a probability density function used to select 

the particle direction, s = | n | is the distance of 

the selected collision point from the previous collision 
point, A(s) is an importance factor for each region which 
changes discontinuous ly at region boundaries, and a(s) is an 
effective cross section which changes discontinuous ly at 
region boundaries and which may change continuously within 
a region. 


* / \ 

The derivative of the logarithm of p^iz^) involves only 
those terms which change when a shield dimension changes, 
i . e . , 




f S 


rco 

ln Pta(2 kn )= £ f 

-1 a(s * )ds • 

- In 

dt | 

1 A( s 1 )a ( s 1 )exp - 


*'o 

o L 


(IT) 


Let s^ denote the distance to a boundary involving the /th 
shield dimension. If the first term on the left side of 
equation 17 is affected by a change in this shield dimension, 
i . e . if s > s^ , then 


at 


i 


- 

•'A 


-I a(s 1 )ds 1 


t \ 

- Ms i ] pjr 


= -»(»,) "iTT 


( 18 ) 


T 7kn 


where a(s^) is the effective cross section at the boundary 
of the shield and is the cosine the particle path makes 

with the outer shield normal. If there is any crossing in- 
volving the j th shield dimension, the second term in equation 
18 will always have a non-zero derivative, i.e., 

9 



St 


In 


i 


CO 


A ( s * )a ( s ' ) exp 


- a 

- 


a(s")ds"ds , 


A(bj, )a(s* ) — — exp f - f 4(s ' )ds 


y' 'r 7 ^ 


^ - J4(s’ )ds' 


j^ A(a ' )a(s ' ) exp t-r a(s")ds"J dsj 


( 19 ) 


Curved shield surfaces may be crossed more than once along 
the path between two particle collision points. Therefore, 
a summation of equations 18 and 19 over every intersection 
involving the ith shield dimension is required to completely 
evaluate equation 17- 


2 . 2 Optimization Procedures 

The shield optimization calculation yields the set of shield 

i * » » . 

dimensions t' = (t^, t 2 , , t L ) such that the dose 

rate, D(t'), meets the dose constraint. The Monte Carlo cal- 
culation is performed for an initial set of shield dimensions 
t = (t^, tg, . . . t^ , . . . , t L ) and yields a set of fluxes, 0^(t), 
3 = 1, 2, ..., J and derivatives, d0j (tK, ^ 2 , ..., J; 

d 1 1 

Jl = 1, 2, . .., L. The assumption is made that the fluxes 
vary exponentially with respect to shield dimension changes 
in the form 


0j(t') = ( t ) 


[ 


exp | a. • ( t 1 



( 20 ) 



where = (a^, a j2 a jL^‘ ^ follows that 


b0At) 



= 0.(t) exp[a. . [a^ . (t'-t)] 


= 0 )^' )a 3£ 


( 21 ) 


In particular 


30,(t) 

at i 



0j(t) 


or 


a 


Jl 


90^(t) / 

3t i / 


0j(t) 


( 22 ) 


(23) 


The weight is also expressed as a function of the shield 
dimensions. The weight is denoted by W(t’) and for spherically 
symmetric shields: 


/ \ 4^ 

W(t') = — 


| p i [ ( r o +t i )3 - r o ] + p 2 [( r o +t i +t 2 )3 -< r o +t i )3t ---} 


Z f,l< r o + Z fc m> 3 - <V Z S m )3 

m=l m=l 


4^ 

3 4 - 7 

£=i 1 


( 24 ) 


il 



where is the density of the j th shield region and r o is 
the minimum shield radius. 


The purpose of the optimization procedure is to minimize the 
weight W(t') subject to the dose rate constraint D(t') = D q 
where D q is a specified dose rate. At this optimum, the 
following equalities hold 



dD (t 1 
9 W ( t 1 


T 


at 


l 


= constant, l = 1, 2, ...,L 


The necessary derivatives are: 




90 (t* ) 

jL 


9 7 


J 


E 


J U 


0j(t) exp 



(t'-t)l(26) 

m 


and for spherically symmetric shield: 



In arriving at the optimum shield, the total shield weight 
is built up in increments of weight AW. Each increment in 
shield weight is always associated with a particular shield 
dimension. At each iteration, the particular shield dimension 
is selected by examining the values of the shield weight quality 
factors, Q^. Each factor represents the approximate change 
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in dose rate per unit change in weight corresponding to a change 
in the ^th shield dimension. Negative Q^'s are the most usual 
and correspond to shields for which an increase in weight -- 
and shield dimensions -- gives a decrease in dose rate. 

Positive Q^'s can occur, however, and correspond to shields 
for which an increase in weight also increases the dose rate. 


If, at a particular iteration, the dose rate is above the dose 
rate constraint, the minimum shield weight increment would 
correspond to the least positive value of those Q^'s for which 
> 0 and for which >t^(min), where t^ (min) is the minimum 
value of the /th shield dimension. If such a exists, the 
dose rate can be decreased while also decreasing the shield 
weight the maximum amount. If there isn't such a Q^, the next 
best procedure is to find the most negative of the Q,^'s for 
which Q^< 0 and for which t£ < t^(max), where t£(max) is the 
maximum value of the j? th shield dimension. A change in that 
would give the maximum decrease in dose rate per unit 
increase in weight. 


If the dose rate is below the specified dose rate at a particular 
iteration, the minimum shield weight increment would correspond 
to the least negative of those Q^'s for which < 0 and for which 
V > ^(roi 11 )* If such a exists, the dose rate can be in- 
creased while decreasing the shield weight the maximum amount. 

If there isn't such a , the next best procedure is to find 
the most positive of those ' s for which > 0 and for which 
tl < t^(max). A change in that §£ would give the maximum 
increase in dose rate per unit increase in weight. 


of the Qi £ 's is selected through 
the above arguments, the corresponding shield dimension t^ 
is changed by a maximum amount At m where At is calculated as 


Assuming a particular value 



AM 

aw (t' ) 

9t i 


(28) 
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I ; • • , , 

If this change, would put t outside one of, its specified limits, 
the value of t would be set to that limit, i.e., t (min)^ 
t m 4:t m (max). , The, shield weight increment AW is calculated as 

D -D(t’) 

o — , 

AW = ( 29 ) 

Q 

m 

subject to the constraint that )AW|<AW q where AW q is a specified 
maximum shield weight increment per iteration. Note that 
AW, and therefore At , may be positive or negative depending 
on the value of and whether the dose rate is above or below 
the dose rate constraint. 

Once a shield dimension is changed, the dose, weight, and their 
derivatives are re-evaluated and the entire process is repeated. 
The optimization would be discontinued in several ways. If 
the dose rate equals the dose rate constraint within the relative 
error of the original Monte Carlo dose rate calculation, the 
program will proceed to the next problem -- which may be identi- 
cal except with more histories to tighten the convergence of 
Monte Carlo calculations. Similarly, if all shield dimensions 
have reached their minimum or maximum values, and if the optimum 
shield cannot be determined with these contraints, the program 
would again proceed to the next problem. Finally, if the dose 
rate and dose rate constraint are decades apart in value, the 
program would reevaluate the fluxes and their derivatives by 
Monte Carlo every time the dose rate changed by more than a 
specified factor during the optimization procedure. 

2 «3 Importance Parameter Optimization 

The optimization of the importance sampling must be performed 
for some function, e.g., dose rate, of the energy-dependent 



fluxes since there is a different optimum for every initial 
particle energy. Therefore, assume that a minimum variance 
calculation of the dose rate is required where 



( 30 ) 


where N is the total number of histories and D is the dose 

_ n 

rate from the nth history and is the average value of the 
dose rate after N histories. The relative error of this dose 
rate is given by 


1 

2 

( 31 ) 



Taking the logarithm of this equation and then performing a 
formal calculation of the partial derivative with respect to 
an unspecified parameter a yields 

aa lnE N = - §~ln D N - In N + ^ In 



_d_ 

d a 




+ 
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Thus the partial derivative of the relative error with respect 
to the parameter a is: 


3 % 

a a 


n2d X 



( 33 ) 


The dose rate from the nth history is given by 




( 34 ) 


where J is the total number of energy groups, k is the number 
of particle collisions, is the flux to dose rate conversion 
factor for the jth energy group, and 0.^ is the flux in the 
jth group from the kth collision of the nth history. Since 


SD N _ 1 y' ar) n 

b "a N a 


( 35 ) 


the calculations required to evaluate equation 33 all involve 
the summation of terms which involve 


k 

The remainder of the analysis, therefore, can be concentrated 
on the partial derivatives of the fluxes. All other operations 
which must be performed are given above. 


a a 
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The fluxes typically depend on the detector position y_, so 
the equation for the particle flux is written as 


^jkn^^ Sjkn^-kn^ 


(3T ) 


The transport kernel K j ( z^ , £ ) does not involve any importance 
sampling parameters so that 


d a 


= ( \n ] K J ( -to’ 2) 


d a 


This equation can also he written as 


a *W 2> 

da 


= s jkn ( Hkr,) 


-kn ' “j v -kn, J 


d 

da 


In 


'WV ) 


(38) 


(39) 


Without going into great detail, it turns out that the particle 
weight S 'i cn (Hi cr) ) is composed of a purely analytical numerator, 
Vjkn^iJkn^ and a denominator which is the product of all the 
probability density functions used to select the collision 
.points, i.e. , 


S jkn^— kn ^ 


V ,1kn^— kn^ 


k 


A V (2 i n) 


Therefore 


In s Jkntekn> 


ln v Jkn ^ ) " ll, /=0 p Jn ( ^n> 


(AO) 


(Al) 


17 



Since v^Cu^) does not explicitly involve any importance 
parameters, it follows 'that 


ln ^ jkn^— kn > = ‘ |r ln p /n<%n> 


da 



d , * / \ 

ST ln V'V 


(42) 


Therefore, equation 39 can be , re-written as 


^.Ikn ( g> 

3a 


k 


^jkn^ S 9a lr) ?j£n^n^ 

J= o 


(43) 


Moreover, the partial derivatives are energy-independent so 
that equation 36 becomes 


3b- Z(Ev^>) (-±£ 

k \H / \ i=0 


ln > ( ^n> 


(44) 


The evaluation of the partial derivatives of the probability 
density functions can be written as 


k 


A ln = 


k-l 


Z^T ln p /r^n )+ ^ ln (“5) 


i=° 
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At the kth collision, the first term on the left side of equation 
45 is known, identically zero if k=0. Therefore, the analysis is 
completed after examining the calculation of the second term. 


At this point it is necessary to identify the particular 
importance parameter a. Since most of the importance sampling 
parameters have fairly involved roles, the technique will be 
applied here to a set of parameters which can have a reasonably 
simple role. These parameters consists of the relative importance 
I of each region. Normally these parameters are all equal. 
However, in asymmetric problems, it turns out that some regions 
are much more important in terms of their scattering contributions 
to a detector. Therefore, these important regions have a 
larger value of I . 

The region importance enters into the selection of a collision 
point through the following probability density function: 


kn 


( ^kr> = 


* 

I r p r 

H 


(s) 


h=l 


Vh 



where r is the region in which the collision occurs (selected 
at random), P r (s) is the piecewise continuous probability 
density function in this region at the selected collision 
point (a distance s from the previous collision point), H is 
the total number of regions in which the collision could have 

"ft -#• 

occurred, and is the integral of P h ( s ') over the partial path 
length in region h. 


Calculating the logarithm of each side of the equation yields: 


ln 


H 

In I r + In p*(s) - In ^I^P* 

h=l 


( 47 ) 


19 



The partial derivative of equation 47 with respect to the 
specific importance parameter . I -- the relative importance 

o 

of region g -- yields 



ln Pkn^kJ 



H , 

E Vgh 

•' h=l 

" H- * 



( 48 ) 


where 6 = 0 if region h is not region g'and 6 = 1. 

gn gg 

Thus equation 48 is evaluated during the random selection of 
the kth collision point and the final term necessary to evaluate 
equation 45 and all preceding equations has been determined. 


The above analysis is used to calculate the partial derivatives 
of the relative error of the dose rate with respect to the 
relative importance I r of each geometric region, and a similar 
analysis is performed for the other importance sampling para- 
meters. The result of the complete Monte Carlo calculation 
is a set of partial derivatives which, for the region importance, 
are given by 


9E 


N 


ai 






(49) 


9 D 

n 

where -g-j — is obtained from equation 44 using equations 45 
and 48. r 


After the calculation is completed, optimal values of the im- 
portance sampling parameters are calculated by requiring that 
the relative error be zero— not actually achieved of course. 
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By a first order expansion 


e n = 0 




<W 


( 50 ) 


where R is the total number of regions. A simple gradient 
analysis says that I - I should be proportional to 
so that 






I 


r 


+ C 




( 51 ) 


where, by substitution into equation 50, 



t 


The program prints the optimum values of I and other importance 
parameters after completing the Monte Carlo flux calculation. 
This analysis is performed for every response function. 

After more experience is obtained with the technique, the 
program could be modified to change these parameters internally 
corresponding to a specified response function. 
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Section 3 

SAMPLE PROBLEM RESULTS 

Two problems were investigated using the shield optimization 
capabilities of the FASTER-III program. Both problems involved 
a spherical reactor-shield configuration and included primary 
neutrons and both primary and secondary photons. 


The two problems were similar except for the power level, 

375 MW and 600 MW respectively. Both problems used a flat 
radial distribution for the primary neutron and photon source 
distribution. The primary photon source included an infinite 
operation equilibrium fission product term. 

The core radii for the two problems were 82.38 and 96.38 cm 
respectively, corresponding to a power density of 4.53 MW/ft . 
Following the core was a 7*62 cm Be reflector; a 5 cm depleted 
uranium shield; three depleted uranium-bora ted water shields 

O 

of 57, 15, and 15 cm thickness and 6.4, 4.6, and 2.8 gm/cnr 
density respectively; and a 117 cm borated water shield. 

This base line shield configuration was based on parameters 
obtained from SANE-SAGE calculations and subsequent calculations 
using the UNAMIT program. Reference 7» The reactor-shield 
compositions are given in Table 1. 

The primary neutron transport calculation utilized multigroup 
cross sections for 26 energy groups. Fifteen energy groups 
were utilized for both primary and secondary photons. The 
secondary production cross sections included both inelastic 
and capture gammas. 

These initial configurations were each analyzed for a point 
detector 30 feet from the core center by following approximately 


23 



TABLE 1 


SPHERICAL REACTOR-SHIELD CONFIGURATION 






COMPOSITIONS (10 24 atoms/cm 3 ) - 


Element 

CORE 

REFLECTOF 

; u 238 
SHIELD 

• 

MIX 1 
SHIELD 

MIX 2 
SHIELD 

_ : - 1 

MIX 3 
SHIELD 

H„0+B 

SHIELD 

H 

0.01976 

0.0 

0.0 

0.0451 

0.0516 

0.0580 

0.0645 

Be 9 

0.0 

0.120 

0.0 

0.0226 

0.0258 

0.0290 

O.O337 

B 

0.0 

0.0 

0.0 

0.000671 

O.OOO 766 

0.000862 

0.000958 

0 

0.01184 

0.0 

0.0 

1 

0.0 

0.0 

0.0 

0.0 

Al 

0.0512 

0.0 

0.0 

0.0 

0.0 \ 

0.0 

0.0 

Zr 

0.01744 

0.0 

0.0 

0.0 

! 

0.0 

0.0 

0.0 

u 235 

0.000979 

- 0.0 

0.0 

0.0 

0.0 

0.0 

0.0- 

U 238 

0.000078 

0.0 

0.0482 

0.01446 

0.00964 

0.00482 

0.0 



500 energy-dependent packets of primary neutrons and photons 
and approximately 7000 packets of secondary photons. The 
dose rates obtained from these calculations are tabulated in 
Table 2 including a breakdown by secondary source region. 

Each of these problems required about 28 minutes on the 
UNIVAC 1108 computer. 

The basic calculated dose rates and dose rate derivatives were 
also used by the PASTER- III program to calculate the minimum 
weight shield configuration which would give a dose rate of 
0.25 mr/hr at the specified detector point. The final shield 
configurations following the optimization are given in Table 3. 

In both cases, the optimum shield configuration is significantly 
different than the base line configuration. Since the base 
line configuration was not generated by the FASTER-III program 
it is difficult to discuss many factors entering into that 
calculation which would account for the different optimal 
configuration. It is noted, however, that the base line 
configuration was generated using parameters corresponding 
to a calculated dose rate an order of magnitude below the 
specified dose rate constraint. Reference 8. As such, the 
base line configuration used in the FASTER-III program was 
determined from an extrapolation of a different base line 
configuration . 

A more critical critique can be made of the FASTER-III results 
independently. First it is noted that neither problem saw a 
significant contribution from photon sources in the core 
region. In fact, the 600 MW reactor dose rate from this source 
was about a factor of two less than it was for the 375 MW 
reactor. This difference is ascribed to the problem statistics 
since core photon sources see approximately 30 mean free paths 
of shield material. Therefore, it is doubtful if this dose rate 
component is converged even with a factor of two after only 
500 packets. 
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TABLE 2 


RESULTS OF BASE LINE CALCULATIONS. 
AT 30 FEET FROM CORE CENTER 



375 MW 

600 MW 

DOSE RATE 

. REACTOR 

REACTOR 

COMPONENT 

(mr/hr) 

Cmr/hr) 


Photon Source Region 
Core 

0.009 

0.004 

Reflector 

3.5xlO -6 

6.3xlO' 6 

Depleted Uranium 

3.2xlO" 5 

1.3xlO~ 5 

Mix 1 Shield 

0.018 

0.026 

Mix 2 Shield 

0.062 

0.075 

Mix 3 Shield 

0.017 

0.063 

Bora ted Water Shield 

0.011 

0.022 


Total Photons 

0.120+0.034 

0.187+0.054 

Neutrons 

0.020+0.002 

0.027+0. 003 

Total 

0.140 

0.214 
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TABLE 3 


RESULTS OF SHIELD OPTIMIZATION 
(O .25 mr/hr at 30 feet) 


Quantity j 

... 

375 MW 
REACTOR 

600 MW 
REACTOR 

Initial 

Fina 1 

Initial 

Final 

Dose Rate (mr/hr) 
Photon 

0.120 

0. 126 

l 

0. 187 

0.153 

Neutron 

0.020 

0.124 

0.027 

0.097 

Total 

0. 140 

0.250 

0.214 

0.250 

Shield Weight (lO^kg) 

I 




Depleted U 

10.2 

12.6 

13.8 

0.0 

Mix I 

71.2 

0.0 

89.2 

6.6 

Mix 2 

22.1 

52.4 

26.4 

52.4 

Mix 3 

16.1 

12.2 

19.0 

63.1 

Water 

86.7 

80.3 

97.7 

85.3 

Total 

206.3 

157.5 

; 246.1 

207.4 

Shield Thickness (cm) 


1 



Depleted U 

5.0 

6.1 

5.0 

0.0 

Mix 1 

57.0 

0.0 

57-0 

7.0 

Mix 2 

15.0 

57.3 

15.0 

48.4 

Mix 3 

15.0 

13.5 

15.0 ' 

5L.4 

Water 

117.0 

120.8 

117.0 

98.4 
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The small contribution from core photon sources decreases the 
amount of high Z shields required around the core. Therefore, 
both problems gave a significant change in the first two 
shield dimensions during the optimization. In the 375 MW pro- 
blem, the first mixture of depleted . uranium-borated water 
(p= 6.4 gm/cm^) was eliminated entirely. In the 600 MW problem, 
the depleted uranium and most of the first mixture were 
eliminated. ■ .. 

The main difference between the two FASTER-III calculations 
was the shift in the .. placement of lighter shield mixes towards 
the core for the 600 MW problem. An examination of the sec- 
ondary photon dose components indicates that the contribution 
from the outer two shields was about 25$ for the 375 MW reactor 
and almost 50$ for the 600 MW reactor. Since these sources 
depend on the neutron attenuation through the closer regions 
and since lower effective Z materials are better neutron 
attenuators on a weight basis, the 600 MW problem tends to 
replace high effective Z material with a lower effective 
Z material. 

The differences in the contribution from secondary sources 
in the outer shield regions is greater than expected for the 
nominal difference in the core region. . Therefore, much of the 
difference in these sources must be ascribed to statistical 
variations. In fact, both problems had approximately 25-30$ 
calculated relative error in tiie total photon dose rate. 

It should be noted that the FASTER-III program includes a 
number of importance sampling techniques which could be used 
to decrease this error. However, both problems were run 
using the built-in definitions of importance parameters. 
Alternatively, more histories could have been used although 
the computer time requirements would have become excessive. 
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Section 4 


CONCLUSIONS AND RECOMMENDATIONS 

The FASTER-III program was developed to calculate neutron and 
photon fluxes at specified points in complex geometries. 
Alternatively, it can also calculate fluxes averaged over specified 
surfaces and volumes. The program was designed such that data 
preparation is simple and so that very little judgment is 
required to set up the importance sampling for most problems. 

The FASTER-III program satisfies these requirements very well. 

The shield weight optimization capability included in the 
FASTER-III program permits the calculation of both base line 
radiation levels and optimal shield thicknesses all in a single 
computer run. However, the very large attenuation factors 
involved in the demonstration problems yielded some questionable 
results. In particular, the statistical differences in the 
relative contribution from various secondary source regions 
caused corresponding variations in the relative distributions 
of shield materials. Of course the statistical variations 
would be less in problems with less overall attenuation. 

The effect of statistical differences on the shield optimization 
can be reduced by following more packets. However, the computer 
times start to get excessive if this is the only approach used. 

It would be more fruitful in terms of the routine application 
of the program to expend some effort towards altering the im- 
portance sampling. 

The FASTER-III program has the capability of calculating optimal 
importance parameters based on partial derivatives of the 
variance. This feature can be used in determining better im- 
portance sampling parameters for shield optimization problems. 

In fact, the overall program efficiency could be improved if 
this feature was utilized on a wide variety of problems with 
the results being used to improve the built-in importance 
sampling models and parameters. 

29 



The relative expense of using the FASTER-III program in a 
somewhat iterative fashion to do the initial sizing of a shield 
configuration should be considered. The least expensive 
procedure suggested for the initial setup of a problem would 
involve either hand calculations and/or point kernel calculations. 
In view of this sizing problem, it is recommended that a point 
kernel option — removal and/or moments data for neutrons, 
buildup factors for photons.,-- be built into the FASTER-III 
program. This option would also include a calculation of 
secondary photon contributions. This option could be used in 
conjunction with the shield optimization procedure and permit 
relatively accurate sizing with an order of magnitude or more 
reduction in computer time when compared with the use of 
Monte Carlo. A very positive advantage of this option is 
that most of the data cards used in a point kernel problem 
would be used directly in the corresponding Monte Carlo problem. 
The incorporation of a point kernel option in the FASTER-III 
program does not involve very major modifications. In fact, 
a similar option was used in one modification of the original 
FASTER program for calculating secondary photon dose rates. 
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Appendix A 
MONTE CARLO METHOD 

The Monte Carlo method as used in the PASTER program is 
described in this appendix. The development starts with the 
order-of-scatter (Neumann series) solution of the transport 
equation. The Monte Carlo method is then applied to the 
spatial integrations. The presentation is of a summary 
nature and no proofs are given. 

1. The Transport Equation 

The particle energy is immediately cast into a multigroup 
framework where the ith energy group includes all particles 
with energies E between the group boundaries E^ and E^ + ^. 

In the conventional manner, higher group indices will indicate 
lower particle energies, i.e., the energy group boundaries are 
monotonically decreasing, E i >E i+ ^. 

The differential angular source of particles in the ith energy 
group which have had exactly k interactions or collisions 
since being emitted from a known independent source is denoted 
S ik ( £'H)* the number of particles in group i per cm^ per 
steradian coming out of a collision at x and proceeding in 
the direction u. 

The differential angular flux of particles in the jth energy 
group due to source particles which have had k collisions is 
denoted by ^ k (£,v), the number of particles in group j per 
cm per steradian crossing a detector at y while heading in the 
direction v. 

The differential angular flux is directly related to the diff- 
erential angular source by a simple line integral over the 
space points which can contribute in the fixed direction v. 
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Allowing for the application of this development to charged 
particle transport, the relationship between flux and source 
is written as: 





CO 




(£-tv,v) K ±J (y-tv,y)dt 



(Al) 


In general, the kernel K^x,;^) is the probability that a 
particle starting at x in group i will arrive at jr in energy 
group j. For neutral particle transport, this kernel is simply: 




exp 



(x+su )ds 



t = 



u = (y-x)/t 


(A2) 


where 2 i (z) is an appropriately averaged total cross section 
for energy group i at the point z. The quantity 5 . . is the 


Kronecker delta function, 
i = j. 


ij 


= 0 if i / 




= 1 if 


The "next collision" angular source, S ± k+1 (x,u), is, in turn, 
determined from the angular flux as: 




j V 


/v s ’- 


) T j ± (x,u.v)d 2 v 


(A3) 


This calculation requires an integration over all initial 
directions v (differential solid angle dfi = d 2 v) which can 
be scattered into the direction u. A summation over all initial 
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groups j for which particles can scatter into group i is also 
required. The kernel Tj ± (x,m) is the probability per unit path 
length per steradian that a particle at x will scatter from 
group j into group i while being deflected through an angle 
cos . 


Monte Carlo Integration 


The Monte Carlo method is used to reduce the integrations 
above to one-point numerical quadratures where the point is 
selected at random. Assuming the most desirable solution 
is represented by the scalar flux at a specified point y, 
then this point solution is composed of contributions from 
the many orders-of -scatter : 


v*> - 

X) v?> 

k 

(A4) 

where <^ k (y) is the scalar flux in group j at y from particles 
which have had exactly k collisions. Each order-of-scatter 
component of the flux can be written as volume integration 
over the kth scattered source since: 

V*> 

* /v 2 '- )d2 - 

(A5) 


= M^/s^Cy-ty^y)^ j ( y-ty,y)dt d 2 v 
4 J 1 o 

(A6) 


~ J f IXkte-^^i j t 2 dtd 2 v 

4k ° 

(A7) 
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2 2 

The differential volume element t dtd v is written in the 

3 ~ ' 

more general form d x where x represents the space point 
£-tv. The directionality of the source is changed from v 
to the more general ir by including an appropriate delta function 
in the integrand. The final form of the equation for the kth 
scattered flux in group j at £ is: 

V £) = S/ S ik ( S.«)K 13 (x, 2 )— a 3 x (AS) 


t = \l-x\ > Z = (z-x)/ t 

This equation is unchanged if the integrand is multiplied and 
divided by the probability density function (pdf) p*(x): 


V £) 


fv 1 

! 

S ik^- , “) 

1 


IldT - 

l 

> 


K ij(£>5L) j P^(x)d 3 x (A9) 


where 


Pv (x) > 0 for all x 


pj(x) > 0 if /s ik (x,u)d 2 u 

1 ** 

yp£(x)d 3 x =1 


> 0 


(A10) 


(All) 


(A 12) 



Equations AID and. A 12 are necessary conditions in defining a 
pdf; it must be non-negative and must integrate to unity. 

The condition in equation All is stronger than is required 
for calculating the flux at % — as written, the same pdf can 
be used for calculating the (k+l)th scattered source for the 
next order scalar flux component. 

The quantity in the first set of brackets in the integrand 
of equation A9 represents a modified source density denoted 
by: 


S ik 


(x,u) 


S i i<-(x*h) 

w 


(A13) 


The second bracketed quantity in equation A9 depends on the 
particular point y at which the scalar flux is being calculated. 

Equation A9 is integrated by selecting a single point z^ at 
random from the probability density function p*(x). The — . 
mechanics of the random selection process will not be discussed 
here. It suffices for this discussion that a value 
obtained at random from pj* (x), gives a one point quadrature 
estimate of the value of the integral in equation A9. This 
estimate is denoted by and is given by 

V 2) = 12 S fk(2k>y) K i;j(Zk>2) ^7?^ (A14) 

i c k 

fc k = l^kl ' *k = (2*lk)/tk 
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Of particular interest is the fact that equation Al4 holds 
for all detector positions y and all energy groups j . It 

• .>r." . ',..■■■* • " . ' ‘ ■ 

is true that there is a particular pdf p£(x). that is best,, 
for calculating the kth scattered scalar flux in group j at 
Zf namely ... .■ ; . 



(A 15) 


Sampling of this pdf to obtain. the discrete point z k will 
give an exact solution for the kth scattered scalar flux in 
group j at y. However, there are several reasons for not 
doing so. Foremost is the fact that in defining the pdf 
through equation A£ it is necessary to essentially calculate 
the flux analytically since the denominator is the desired 
answer. Second, it is virtually impossible to define and then 
sample a pdf as complicated as equation AI 5 . Furthermore, 
it is usually more economical to approximate equation AI 5 
for the dominant energy group and then use the selected 
point z k in the calculation of fluxes for all the energy groups 
simultaneously as implied in equation Al4. Finally, this 
optimal pdf only holds for the k th scattered flux, and the 
next step after calculating the kth scattered flux is to 
use the source strengths S ik (z k ,u) in defining the (k+l)th 
scattered source for the next order of scatter. There is a 
different and much more complicated prescription for defining 
an optimal pdf to be used in the selection of for this • 
next step. 


After these negative comments, there are several features 
of equation A15 which are quite useful. It does provide insight 



into the functional form of the optimum Last collision pro- 
bability density function. Furthermore , the more complex pdf 
to be used for future orders-of -scatter can be approximated 
by equation A 15 by simple alterations of the transport kernel 

p 

K ij (*■;£)• Finally, it does include the l/t singularity in 
a manner which obviates any difficulties in calculating 
accurate, finite variance fluxes at a point. 

As discussed above, the pdf pj*(x) is used to select a discrete 
point which is then used via equation Al4 to estimate scalar 
fluxes at the point The point source defined at can 
actually be used in the estimate of angular fluxes since 
the source strengths k J — ^ can be evaluated for various 

directions .Thus, the angular flux at any point can be 
estimated, using these same sources, as: 


•j k (z,v) = ^SJ k (z k ,v) (A 16) 


s k * |z-5k| > - (rVAic 


Of course, if the point jr lies in the k th scattered source 
volume, a pdf which includes a dependence on £ should be used 
in selecting z^. However, if y lies outside this source volume, 
there is no difficulty. 

In addition, these point value flux estimates can be area or 
volume averaged. The equations are simply 
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V A <: 


il = iiggd. 

i .a t. 


dA (A 17) 


and * ilc (V,v) = E V ^,(^,2) aV 


E V ik'-=k’- 
1 V 


ij -k J 


(A18) 


In both cases, the integrations are transformed to spherical 
coordinates so that dA = t\ d 2 v k where „ is the normal to 

IV 2 1 


2 2 

the area, and dV = t^ d v-^dt^ so that 


V A 


* Ek f s Ik<5k’i> K ij'V2k +t k ( ^> 

1 V I— k* -I 


= Si s ik<2k>2) K i.1 tl k 


(A 19) 


2k- a 


■(v) 


* ik (V,v) = 


][v f I S ik ( V- )K i j ( ^k^k +t k^ ) fc kd 2 v k dt k 

1 4^t k (v) tk 
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■ , ,1' i. > . 

• •’ % 

: '■ 1 i .* ... 1 " 


E 


v s L ( 2k-^ ) 


(A 20) 


> ''> >v I . 


where t k (v) is the distance from along v. up to the area 
or volume, and - ^(v) is the .distance through the '.volume..! 

Any summations over /mu ltiple the line z,+t, v 

with the area or volume are implicit. In; 'the' above equations ^ ’ ■ 

The point angular f luxes defined by equation Al6 are used to 
define the (k+l)th scattered source: ■ , , [ , : 


S 1)k+1 (5,H) = Tjifx.v.n) :i r v 

j , r;' : 

This reduces to: * jV— 


(A21) 


S i,k + 1 ( *^ , = .^jk^^k^^i^VH)/ 

.*. ■■ .. ' j >. :/ .i . ' 


(A22) 


since the angular fluxes are defined ’aC^feifag" •in' , --;the.>dire 5 .ction : . 
v k only, where v k = (x - z k )/.| x - z k | . / ■ ^V. 


As before, a pdf y p^ + ^(x)., the (k+l)th source is.' defined and 
then sampled to obtain... a discrete value z k+k of x. The source 
strength at this point is then denoted by s * k+p^-k+l’-) and is 
given by : " ■ : ‘ ■ 1- ' 
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(A23) 



k+l^+l'-^ 


S i ,k+I^-k+L J 
” p k+l^k + l 


H) 

T 


~ y^^1k^ g k+l ,Y k^ T ji^~ k+l^k* 11 ^ 
J Pk+l^k+J 


(A24) 


Since flux estimates may require evaluation of S.^ k+L (z k+1 ,u; 
for various directions u, it is expedient to define point mono- 
directional values of the flux going into the point z k+1 : 


*Jk (Sk+l'V = < ^ik^-k+l ,Y k ) 

Pk+l^k+l^ 


(A25) 


The angular point source S* Q is determined from input 

only for the independent source. In all other instances, it is 
determined from the equation 


S i,k+l^-k+l , ~^ " ^k^k+l ,Y k^ T i^-k+l ,Y k ,:u ^ ( A26 ) 


The process of reducing each volume -distributed order-of -scatter 
source to a poibt representation by random sampling, of using 
the same point representation to give the volume-distributed source 
for the next order of scatter is continued until all the 
particles at a given order of scatter no longer yield a 
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significant contribution to the flux estimates. The estimate 
of the total angular flux at a point y from this procedure 
is then obtained by summing the individual order-of -scatter 
components : 




2J V 2 '-’ 

k 

23 £ S !k ( -k’ L> 
k i t 2 


(A27) 


(A28) 


The total scalar flux is obtained by a simple integration over 
solid angle which yields 


Vz) = £ K i;l f 2 k ^ ) < A2 9) 

k i .2 

r k 

This process gives a single, inaccurate, estimate of the total 
flux. Therefore the process is repeated a specified number of 
times and the average of all the estimates is accepted as the 
best estimate of the total flux. Introducing the subscript 
n to denote the iteration index, then <l>. (y) is the estimate 
of the flux in the jth energy group at y of particles obtained 
on the n th iteration. If N is the total number of iterations, 
then the total flux is estimated by: 

N 

Sjte) = ST (A30) 

n-1 
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Since each iterant represents an independent estimate of the 
flux, it is possible to approximate the standard deviation 
of the total flux by: 
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(A31) 



